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We report a first-principles theoretical study of hyperfine interactions, zero-point effects and defect energetics of 
muonium and hydrogen impurities in silicon and germanium. The spin-polarized density functional method is 
used, with the crystalline orbitals expanded in all-electron Gaussian basis sets. The behaviour of hydrogen and 
muonium impurities at both the tetrahedral and bond-centred sites is investigated within a supercell approxima- 
tion. To describe the zero-point motion of the impurities, a double adiabatic approximation is employed in which 
the electron, muon/proton and host lattice degrees of freedom are decoupled. Within this approximation the 
relaxation of the atoms of the host lattice may differ for the muon and proton, although in practice the difference 
is found to be slight. With the inclusion of zero-point motion the tetrahedral site is energetically preferred over 
the bond-centred site in both silicon and germanium. The hyperfine and superhyperfine parameters, calculated 
as averages over the motion of the muon, agree reasonably well with the available data from muon spin resonance 
experiments. 
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I. INTRODUCTION 

Hydrogen is known to have a wide range of physi- 
cal effects in semiconductors, including the passivation 
of states associated with deep-level impurities, enhance- 
ment of the diffusivity of oxygen, and the formation 
of large, planar structures known as platelets. It is 
present in large quantities during the processing stages 
of device manufacture and is one of the commonest im- 
purities in technologically important materials such as 
silicon and germanium. Since such hydrogen impurities 
can have significant effects on semiconductor electrical 
properties, a more complete understanding of their be- 
haviour at the microscopic level is desirable. 

Paramagnetic hydrogen centres can in principle be 
studied using the electron paramagnetic resonance 
(EPR) technique. Information on their local environment 
is obtained by following the time evolution of the signal 
corresponding to the coupling of the spin of the impu- 
rity with an external electromagnetic field. However, few 
studies have been reported for hydrogen in semiconduc- 
tors because the hydrogen atoms are mobile and diffuse to 
defects where they form passivated complexes. The tran- 
sient centres of isolated hydrogen impurities are neverthe- 
less of significant interest because of their involvement in 
diffusion processes, and in fact they may be studied using 
muon spin resonance (/iSR) techniques. Muons have the 
same charge as protons but only about one ninth of the 
mass. They can capture an electron to form a hydrogen- 
like bound state known as muonium (given the symbol 
Mu), and it is thus possible to consider the muon as a 
proton analogue. Transient centres of implanted positive 
muons in semiconductors may be studied as the muon 
has a lifetime of just 2.2 /is and diffuses to locally stable 
sites within a few nanoseconds. The short lifetime also 



means that there is almost never more than one muon 
in the sample at any one time, and that the distribution 
of muons does not reach true thermal equilibrium. In a 
muon spin resonance experiment fully polarized positive 
muons are injected into the sample, and by observing the 
positrons produced by the decay of the muons one can 
obtain information about the defect. 0] 

When /iSR experiments are performed on silicon or 
germanium, two different hyperfine signals are observed. 
One of these is entirely isotropic while the other has an 
anisotropic (dipolar) component with uniaxial symmetry 
along the [1 1 1] axis. [Q The impurity responsible for the 
former signal is usually referred to as normal muonium 
(Mu) and that for the latter, anomalous muonium (Mu*). 
Normal muonium has been identified as muonium in the 
interstitial region, probably in the vicinity of the tetra- 
hedral (T) interstitial site. Symons and Cox j2j first sug- 
gested that anomolous muonium corresponds to a neu- 
tral muonium at the bond-centred (BC) site and this has 
been borne out by a number of theoretical studies. The 
various experimental data for muonium in silicon have 
been interpreted in terms of a configuration-coordinate 
diagram. ^ 

The majority of recent theoretical work in this area 
has been at the first principles level within an adiabatic 
approximation, using the local spin density (LSDA) or 
generalized gradient (GGA) approximations to density 
functional theory (DFT). B Calculations using pseu- 
dopotentials and plane-wave basis sets 1^-0] have been 
reasonably successful in reproducing the hyperfine and 
superhyperfine parameters observed in experiments. The 
majority of such calculations appear to demonstrate that 
hydrogen impurities at the T and BC sites have similar 
energies. 

The application of the Feynman path-integral H 
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method to the study of these systems enables the effect 
of the quantum nature of the muon to be studied directly 
at finite temperatures. However, the large computational 
demands of such an approach have limited its use to date. 
Ramirez and Herrero used the path integral molec- 
ular dynamics method to study hydrogen and muonium 
in silicon with the H/Mu-Si interaction described by an 
empirical three-body potential. However, the results ap- 
pear to be in conflict with experiment. Recently, Miyake 
et al. applied the path-integral Monte Carlo tech- 
nique to the study of hydrogen and muonium at the T 
site in silicon, with the electron-electron interactions de- 
scribed within the LDA. Despite finding the T site to be 
a local maximum on the potential energy surface, they 
found the muon distribution to be peaked at that site 
because of the quantum motion. 

In this work we employ all-electron DFT calculations 
within a double adiabatic approximation to study muo- 
nium and hydrogen at the BC and T sites in silicon and 
germanium. The use of all-electron calculations allows 
an assessment of the accuracy of the correction proce- 
dures which are used to obtain the hyperfine and super- 
hyperfine parameters in pseudopotential calculations. 
The use of a double adiabatic approximation allows us to 
obtain both the zero-point energy and wave function of 
the impurity. Our inclusion of the zero-point motion is 
at a level beyond that in previous first principles calcula- 
tions since the positions of the host silicon or germanium 
atoms are allowed to relax in the presence of the zero- 
point motion of the impurity. At this level of approx- 
imation the relaxations of the host lattice are different 
for a muon and a proton. Our calculations thus allow an 
assessment of the differences in the potentials felt by the 
two impurities, thereby testing one of the assumptions 
underlying the configuration-coordinate diagram [| used 
to interpret experimental data. 



II. METHOD 

A. All— electron spin— polarized LSDA— DFT 
calculations 

All the first principles calculations reported here were 
performed with the CRYSTAL95 software package p2| . 
The (zero temperature) spin-polarized density func- 
tional method was used, together with both lo- 
cal density and gradient corrected approximations to the 
exchange-correlation functional (namely the Perdew- 
Zunger LSDA [|5) and the PW91 form of the GGA 
The calculations were performed within a periodic su- 
percell approach with a single hydrogen impurity in su- 
percells containing either sixteen or fifty-four silicon or 
germanium atoms. Fig. |l| shows the relaxed atomic en- 
vironments of a single muon at both bond-centred and 
tetrahedral impurity sites in silicon. The measured lat- 



tice constants (5.429 A for silicon and 5.6579 A for ger- 
manium) were used in all calculations. 

Other approximations made were as follows. The use 
of local basis functions requires the real space Coulomb 
and exchange series to be limited and approximated as 
described in references jljjl^; the accuracy with which 
the various Gaussian integrals are computed is controlled 
by classifying basis function pairs according to overlap 
or penetration criteria defined by five parameters, which 
in this study were set to IQ-"^, 10"*^, 10""^, 10""^ and 
10^^'*. This is normally sufficient to give a numerical 
error of less than 0.001 eV/atom in the relative ener- 
gies of different structures. The reciprocal space integra- 
tions necessary to reconstruct the density matrix in real 
space at each self-consistent cycle were approximated by 
summing over a set of k-points belonging to a mesh of 
Monkhorst-Pack ||l8| type which was centred on the ori- 
gin in reciprocal space. The convergence of both the 
total energy and the isotropic hyperfine parameter of the 
muon with respect to the reciprocal space sampling den- 
sity was investigated. A 4x4x4 k-point mesh was found to 
be sufhcient for the 16-atom supercell. With this mesh, 
the total energy and isotropic hyperfine parameter are 
within 0.0025 eV/atom and 3 MHz of their fully con- 
verged values, respectively. A 3x3x3 k-point mesh was 
used for the 54-atom supercell, which also gives excellent 
convergence. The convergence of various quantities with 
respect to the supercell size and basis set is discussed in 



Section III 



A hydrogen impurity introduces a defect state into the 
band gap of the host crystal. Finite supercell sizes give 
rise to interactions between the defects in different cells 
and thus to a small but potentially significant disper- 
sion in the defect band. This dispersion could lead, for 
example, to overlap of the majority spin defect band 
with the minority spin defect band and/or the silicon 
valence/conduction bands. In either case an unphysical 
conducting state is formed. This problem is not entirely 
eliminated even with the use of the larger 54-atom su- 
percell. However, judicious use of the level-shifting con- 
vergence technique [|l^,|l9| allows a small decoupling of 
unoccupied and occupied states which prevents the sys- 
tem entering a conducting state. Population analysis of 
the final self-consistent wave function revealed that each 
supercell contained a single extra majority spin electron 
as expected on physical grounds. The calculations thus 
correctly model this aspect of the behaviour of a single 
impurity in a large crystal, which is necessary in order to 
obtain physical hyperfine and superhypcrfinc parameters. 



B. Gaussian Basis Sets 

The Bloch functions required to expand the Kohn- 
Sham orbitals in the solid-state band structure problem 
are built from periodic arrays of atom-centred Gaussian 
functions. One motivation for the use of such a basis set 
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is that all electrons in the system may be treated explic- 
itly, allowing the spin density at and around the nucleus 
(and hence the hyperfine parameters) to be calculated di- 
rectly without resorting to correction procedures of the 
sort required in pscudopotential calculations. 

The basis set used for the majority of the sili- 
con calculations was of the type s{8)sp{8)sp{3)sp{l) 
where the numbers in brackets refer to the number of 
contracted primitive Gaussians making up each shell. 
For convergence checking we also used a higher qual- 
ity silicon set with an additional d polarization func- 
tion of the type s{8)sp{8)sp{l)sp{l)sp{l)sp{l)d{l). The 
basis set used for the germanium calculations was 
s(9)sp(7)sp(6)sp(3)sp(l)d(6)d(l). 

To describe the hydrogen atom, an uncontracted ba- 
sis of eleven s functions and a single p function was used. 
Such a large set (mainly consisting of functions with rela- 
tively high exponents) was found to be necessary to com- 
pute accurate hyperfine parameters. A spin-unrestricted 
Hartree-Fock calculation of the total energy of the free 
atom with this basis gave —0.49988 Ha which is close to 
the exact result of —0.5 Ha. The isotropic hyperfine pa- 
rameter was 1421.9 MHz compared with that obtained 
from the exact wave function of 1422.8 MHz. The cor- 
responding values obtained from an LSDA-DFT calcu- 
lation with this basis set were —0.47833 Ha (which is 
very close to the value of —0.47885 Ha obtained from an 
atomic code using integration on a very fine grid) and 
1356.6 MHz. 

Optimal Gaussian basis sets for use in close packed 
solids are significantly different from those appropriate 
to the atomic and molecular cases. In particular, careful 
optimization is required to avoid the problems of linear 
dependence and basis set superposition error due to the 
overlap of diffuse functions. In this study we used the 
following procedure. All basis set parameters were first 
optimized in the free atom. The exponents and contrac- 
tion coefficients of the valence functions in silicon and 
germanium were then reoptimized in the pure bulk ma- 
terial. Finally, a hydrogen atom was inserted at a bond- 
centred site, the positions of the nearest-neighbour sil- 
icon/germanium atoms relaxed, and the parameters of 
the valence functions of each atom again reoptimized. 
To test the transferability of the optimized basis sets the 
hydrogen was displaced from the BC site along the bond 
by 0.27 A, and the basis function parameters were reop- 
timized for the new geometry. The energy as a function 
of displacement along the bond was calculated for each 
of these two basis sets. The variation in energy was es- 
sentially the same. 

The final exponents and contraction coefficients of all 
the basis sets employed in this study are available else- 
where. @ 

It is important to investigate the possibility of basis 
set superposition error (ESSE) in defect energetics cal- 
culations for a system described by a localized basis. The 
basis sets for the host lattice atoms are necessarily incom- 
plete. Insertion of an impurity atom allows additional 



variational freedom in the description of the atoms ad- 
jacent to the defect site, particularly when the impurity 
is described by a relatively diffuse basis set. This can 
distort the relative stabilities of defects at impurity sites 
of differing local coordination number and geometry. In 
the present case, the hydrogen impurity is considerably 
closer to its neighbours at the BC site than at the T site, 
and thus one might expect the BC site to be artificially 
stabilized with respect to the T. 

This expectation is confirmed by an estimate of the 
BSSE using a counterpoise correction. For the 16- 
atom supercell, addition of "ghost" hydrogen basis func- 
tions into the relaxed silicon lattice lowered the energy 
per cell by 0.199 eV (BC site) and 0.068 eV (T site) with 
the smaller silicon basis, and by 0.0533 eV (BC site) and 
0.0243 eV (T site) with the larger silicon basis. In ger- 
manium, the energy is lowered by 0.251 eV (BC site) and 
0.0534 eV (T site) by the same procedure. Inclusion of a 
lattice of "ghost" silicon/germanium functions around a 
hydrogen atom lowered the energy by less than 0.0005 eV. 
These numbers may be taken to give a rough indication 
of basis set incompleteness in each case. It can thus be 
concluded that in silicon the BSSE lowers the energy of 
the BC site over that of the T site by around 0.13 eV with 
the smaller basis, but by only 0.03 eV with the larger set. 
The corresponding correction for the germanium case is 
0.20 eV. 



C. Calculation of zero— point motion 

For the calculation of the zero-point motion of the 
muon/proton a double adiabatic approximation was 
used. This means that the motions of the electrons and 
of the muon are considered to be decoupled from the 
motion of the atomic nuclei, and furthermore that the 
electronic motion is decoupled from that of the muon. 
The approximation is justified by the fact that a muon 
is roughly 207 times more massive than an electron and 
around 243 times less massive than a silicon nucleus. For 
a proton the equivalent factors are respectively 1836 and 
28; the decoupling of the proton and silicon motion is 
thus somewhat less justified. The mass differences are 
of course more favourable for the heavier germanium nu- 
cleus. 

The positions of the silicon/germanium nuclei are de- 
noted by r„, the muon or proton positions by r^, and the 
electron positions by rg. The double adiabatic approxi- 
mation is used to decouple the motions of the particles by 
approximating the wave function as a product of nuclear, 
muon/proton and electronic parts, 

*(re,r^,r„) = ^A"(r„)X''(r^;r„),^^(re;r,„r„) , (2.1) 

where the variables to the right of the semi-colons appear 
as parameters and those to the left are dynamical vari- 
ables. Within the double adiabatic approximation the 
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three wave functions each satisfy separate Schrodinger 
equations: 

i?e(re;r,,,r„)0''(re;r^,r„) = ^''(r^, r„) (/)''(re; r^, r„) , 

(2.2) 



i7^(r^;r„)X4'(r^;r„) = i?^i(r„)X^(r^;r„) , (2.3) 



iI„C(rn)-KC(r„) 



(2.4) 



We then performed a further twelve LSDA calculations 
as a function of the position of the muon/proton for each 
of these nuclear configurations in order to map out the 
required potential energy surfaces. A similar procedure 
was used for germanium. For the T site the static re- 
laxation of the NN atoms is very small and we assumed 
that the zero-point motion of the muon/proton would 
not give any additional relaxation. 

In order to solve Eq. 2.3 for the muon/proton wave 
function we used a fitted polynomial for the energy 
i?^(r^,r„). For the BC site we use a cylindrical coor- 
dinate system with the origin at the BC site and the z- 
and /9-coordinates directed along the bond and in the 
plane perpendicular to the bond, respectively. We have 
neglected the Q dependence of the potential. This as- 
sumption was checked by displacing the muon by 0.53 A 
from the BC site along the [-1 1 0] direction and then ro- 
tating it about the [1 1 1] axis. The maximum variation 
seen in the energy of the 16-atom supercell during the 
rotation was just 0.002 eV. The Taylor expansion of the 
cylindrically symmetric potential, neglecting sixth-order 
terms and higher, is 

Vbc(p, z) = yBc(0, 0) + 13 + -iz" + ip'z' + Cp^ -t- T]z^ . 

(2.8) 

As a further simplification in order to avoid costly, low- 
symmetry calculations, the bp^z^ term was neglected. 
The resulting Schrodinger equation is separable. In order 
to check the assumption of p-z separability, a few calcula- 
tions were performed with the muon at points where both 
the p and z coordinates were non-zero. These energies 
where T is the kinetic operator and Kb is the Ewald inter- ^^^^ ^j^^j^ compared with those predicted by the fitted 
action between particles of types a and h. The term V^,^ is potential neglecting the term in p'^z^. The errors due to 



The subscript a labels the different eigenstates of the 
muon. Although only the ground state of the nuclear 
wave function is considered here, it is also labelled by 
a since it depends on the chosen muon eigenstate. The 
different electronic eigenstates are not labelled, since it 
is only the ground state of the electronic wave function 
as a function of the muon and nuclear positions that is 
of interest in the current work. The three Hamiltonians 



are: 



i?e(re;r^,r„) = T,(y,) + 14e(re) (2.5) 
Ve„(re,r„) -f Vep(re,r^) , 

^^(r^;r„) = f^(r^) + F^^(r^) + i;,„(r^,r„) (2.6) 
+ £;^(r^,r„) , 



i?„(r„) = K„(r„)+^^^(i-„) 



(2.7) 



a constant describing the interactions between the impu 
rity atoms in the different supercells. The electronic en- 
ergy, i?*^, appears as a n effective potential in t he muonic 
Hamiltonian, Eq. 2.7. Hence, when Eq. 2.3 is solved. 



the resulting energy includes the electronic contribution. 
This energy then appears in the nuclear Hamiltonian as 
an effective potential. Thus the total energy of the sys- 



tem is given by the eigenvalue in Eq. 2.4 



neglecting the (P'z^ term increased only slowly away from 
the BC site, and the corresponding error in the ground 
state energy is small because the wave function of the 
muon/proton is localized around the BC site. 

Our approximation formula was then obtained by a 
least-squares fit of the twelve energies calculated at dif- 
ferent values of (for fixed r„) to the resulting polyno- 



In order to find a good starting point for the BC 
calculations we performed LSDA calculations with the 
muon/proton fixed at the bond centre. The positions of 
the nearest and next-nearest neighbour (NN and NNN) 
silicon/germanium atoms were then relaxed. For the T 
site the muon/proton was held fixed at the T site while 
the four NN silicon/germanium atoms were relaxed in 
the radial direction. 

The next step is to calculate the potential experienced 
by the muon/proton in the c ryst al by solving the elec- 
tronic Schrodinger equation ( 2.2) within the LSDA, as 
a function of the parameters and r„. For the BC site 
calculations, the parameters r„ were varied by consider- 
ing four different additional relaxations of the NN silicon 
atoms along the [1 1 1] direction. For each of the positions 
of the NN silicon atoms, the NNN atoms were relaxed. 



mial. The fitted values of the parameters in Eq. 2.8 used 



in the calculations of the zero-point energies for silicon 
and germanium are given in Table ^. 

The polynomial expansion for the energy surface 
around the T site is the Taylor expansion which is in- 
variant under all of the 24 rotations forming the group 
Td- Cartesian coordinates centred at the T site were used, 
and in order to obtain a good fit to the energy surface it 
was found necessary to include terms up to sixth order: 

Ft (a;, y, z) = yT(0, 0, 0) + a2[x^ + + z^) + as ^yz 

+ a^{x^ -I- / + z^) + W[x^y^ + x^z^ + z^) 
+ a5 xyz{x^ + + z^) + ae{x^ + / + z^) 
+ beix^y* + x^z* + x^y^ + x*z^ + y^z^ 
+y^z^) + cex^y^z^ . (2.9) 
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The fitted values of the parameters in this equation are 
given in Table ||. 

Th e solution of the muon/proton Schrodinger equation 
( 2^) was performed by diagonalizing within a basis of 
harmonic oscillator eigenfunctions centred on either the 
BC or T site. We constructed muon and proton basis 
sets consisting of the solutions of the harmonic part of 
the calculated potentials: 



and 



K)Bc = ^Bc(0, 0, 0) + p{x^ + y2) + ^ (2.10) 



14)t = V^T(0,0,0) + a2(x2+y' + z2) . (2.11) 



For a harmonic potential, the Schrodinger equation may 
be solved by separation of variables: 

X{v^) = X{x)y{y)Z{z) . 

This gives three separate equations for X, and Z, 
whose solutions are of the form: 



where A, 



Xi{x') = AiHi{x')e--2-'\ (2.12) 
: ^ is the normalisation factor, x' = 

V2'7r2;! 

{2MCx)^x is a rescaled variable which allows the eigen- 
functions to be written in terms of the standard Hermite 
polynomials, Hi{x'), and is the appropriate harmonic 
coefficient. The associated energy eigenvalues are 



1 + 



1 



2Cx 
M 



(2.13) 



where A runs over the three directions, x, y, and z, and 
C\ is the harmonic coefficient corresponding to that di- 
rection. The value of the mass, M , of the particle de- 
pends on whether we are solving for the proton or muon 
wave function. 

Having constructed our basis functions we now con- 
sider the full potentials which are written as the sum of 
harmonic and anharmonic terms: 



^BC = Vbec + AT/bc 
Vt = Vo^ + AVt . 



(2.14) 
(2.15) 



The Hamiltonian matrix elements were calculated in the 
basis of the harmonic solutions and the resulting ma- 
trix equations were diagonalized. A basis set constructed 
from all Hermite polynomials up to and including eighth 
order and containing a total of 729 basis functions was 
found to be sufficient to obtain converged values for at 
least the lowest six eigenvalues, i?^(r„), of the system 
with the muon/proton at the BC site. At the T site it 
was found necessary to increase the size of the basis set 
to include all Hermite polynomials up to twelfth order, 
which gave 2197 functi ons. 

The solution of Eq. 2.7 is trivial as the operator is 



multiplicative and the eigenfunctions are delta functions. 
The total energy, E^, is therefore the sum of E^{rn) and 
the Ewald energy of the lattice of host atoms, Vnn{j^n)- 



D. Hyperfine and superhyperfine parameters and 
motion averaging 

The components of the hyperfine tensor. A, define the 
spin Hamiltonian for the hyperfine interaction between 
the spins of an electron and a nucleus: 



(2.16) 



The hyperfine tensor is normally split into isotropic and 
anisotropic parts. 



A = + A, 



(2.17) 



where I is the (3 x 3) unit matrix. The isotropic hyperfine 
parameter (or superhyperfine parameter when it is calcu- 
lated at one of the nearest neighbours of the impurity) is 
given by 



= 104.9827>^(r„) [MHz] 



(2.18) 



where fiQ is the permeability of free space, /x"^ is the Bohr 
magneton, /i" is the nuclear magneton and ge and gn are 
the electron and nuclear g factors. The position of the 
nucleus is denoted by r„ and = — [bohr~^]. |2|] 
The anisotropic part of the hyperfine tensor is given 

by 

= £geM'5nM" J T(r)p,(r + r„)dr (2.19) 
= 12. 531 J T{r)p„{r + rn)dr [MHz] 



where T(r) is a traceless tensor. 



1 



3x^ — r^ 3xy 



3xz 



T(r) = ^ 3xy 3y' - H 3yz , (2.20) 



3xz 



3yz 3z — r 




and the origin of coordinates is on the nucleus at r„ . For 
a particle located precisely at the BC site the hyperfine 
interaction has axial symmetry with respect to the [1 1 
1] axis and thus Ap has the form 



(2.21) 



with Ap being the anisotropic hyperfine parameter at this 
site. For a particle located precisely at the T site, all 
elements of Ap are zero and hence the hyperfine tensor 
is purely isotropic. 

In reality the muon/proton will explore the environ- 
ment around these sites by virtue of its zero-point motion 
and thermal effects. In order to account for the zero- 
point motion the hyperfine interaction tensor must be 
averaged over the squared modulus of the muon/proton 
wave function: 
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(A), = 



J |X(r^; r„)|2A(r^)dr^ . (2.22) (A^p)^ = J^Y. j j P^^^) O^C^m; r«) 



To evaluate the integral for each component of A we fit 
As{x^, Hfj^jZfj) and each of the six distinct elements of the 
symmetric tensor Ap{xp,,yp,, z^) to polynomial expres- 
sions of the correct symmetry. Since the muon/proton 
wave function is expanded in terms of Hermite polyno- 
mials, analytic expressions for the elements of (A)^ may 
be obtained. 

For the isotropic hyperfine parameters the polynomial 
expression for As(x^, y^, z^) has the same symmetry as 
the relevant potential energy surface. These parameters 
were expanded in sixth-order polynomials. The polyno- 
mial describing the isotropic superhyperfine parameter 
at the BC site contains terms which are odd in z^. (Su- 
perhyperfine parameters were not calculated for the T 
site.) Each of the elements of A J and A^^ were fitted 
to second-order polynomials of the correct symmetry. 

We now consider the symmetry of the hyperfine ten- 
sor including the effects of zero-point motion. Within 
the double adiabatic approximation the muon motion is 
described by the wave function X(r^;r„). The motion- 
average of the a/3-component (a,/? = x,y,z) of the 
anisotropic hyperfine tensor is given by: 

{Aai3)t,^C j j |X(r^;r„)|Va(r-hr^)T„^(r)drdr^ , 

(2.23) 

where C is a constant, is the electron spin density, 
and T gQ de notes the components of the tensor T defined 



in Eq. 2.20 



The muon/proton may be said to be trapped in a po- 
tential well if its wave function is negligibly small out- 
side of an equi-potential-energy surface enclosing the re- 
gion. The muon wave function, X, is the non-degenerate 
ground state of the potential well and therefore has the 
full point-group symmetry of the well, i.e., 

P{Qi)X{Y^-Yn) = XiH-^r^^-v,,) ^ X{Yf,-Y^) , (2.24) 

where P is a scalar transformation operator, Qi is an 
operation of the point group of the well, and is the 
corresponding transformation matrix. The electron spin 
density, pcr(r + r^), satisfies 

P{Q^)Pa{v + r^) = Pa{^i\r + Y ^)) = (r + r^) . 

(2.25) 

{Aap)^ is unchanged by a scalar transformation of the 
integrand, i.e.. 



(2.26) 



{A<,p)>.^C j j P{Qi)[\X{Y^-Y,,)\'x 

Pa{Y + Yf,)Tai3{Y)] d-Y d-Y 



Pa{Y + r^)TQ^(r)] fir fir p 



(2.27) 



Using Eqs. 2.24 and 2.25 we have 



(^a/3)p--§/ / |X(r^;r„)|^p,(r + r^) X 



N 



dY dY,, 



(2.28) 



3erties of (^a/s)/^ are easily ob- 
For example, {Axy)fj, will be 



{Aai3)fj, is again unaltered if we sum over the i operations 
and divide by their number, iV, 



The symmetry pro 
tained from Eq. 2.28 
equal to (Ay,)^ if J2,P{Qr)xy = E^P{Q^)yz■ Tak- 
ing the specific case of the T site, it is easily shown that 
J2i P{Qi) = J2i P[Qi) = where the sum is over 
the 24 operations of the tetrahedral point group. Similar 
arguments show that all the elements of {Aap)^ are zero 
for the T site. Similarly, for the BC site we find that 
all off-diagonal elements of {Aap)^ are equal, and the 
diagonal elements are zero. 

If the zero-point motion of the muon is neglected then 
|X(r^; r„)p = 5(r^ — Tq), where ro is the position of the 
muon. It follows that if the muon is placed at an invari- 
ant point of the symmetry group of the well, then includ- 
ing the zero-point motion does not change the symmetry 
of the anisotropic hyperfine tensor. This result explains 
why the zero-point motion does not affect the symmetry 
of the anisotropic hyperfine tensor for either the T or BC 
sites considered here. 

The presence of the muon could lead to a symmetry 
lowering distortion of the host lattice, in which case the 
appropriate point group is the lower symmetry one. We 
have not considered the possibility of symmetry lowering 
distortions in our calculations because of the computa- 
tional cost of evaluating the energy E'^{y^, r„) of Eq. 2.2 
for the required atomic configurations. However, we be- 
lieve such distortions to be unlikely for the cases consid- 
ered here. 



III. RESULTS 

A. Static relaxations 

The static relaxations (neglecting zero-point motion) 
are, of course, identical for the muon and proton. Cal- 
culations with the muon/proton fixed at the BC site 
of the 16-atom silicon cell showed that the two near- 
est neighbours of the muon/proton relax outwards from 
the muon/proton by 0.40 A along the [I I I] axis with 
the NNNs relaxing by 0.01 A in the same direction. The 
corresponding relaxations for the 54-atom supercell were 
0.39 A and 0.02 A. These values are close to the plane- 
wave pseudopotential results of Luchsinger et al. who 
obtained relaxations of 0.45 A for the NN silicon atoms 
and 0.07 A for the NNNs. For the 16-atom germanium 



6 



cell the corresponding relaxations were 0.44 A for the 
NNs and 0.02 A for the NNNs. Our NN relaxation is in 
good agreement with the value of 0.42 A calculated by 
Vogel et al. || 

At the T site, the NN atoms in the host lattice were 
allowed to relax in the radial direction. The relaxations 
in silicon and germanium were approximately equal and 
very small; just 0.02 A towards the muon/proton in the 
16-atom supercell and 0.03 A (in the same direction) in 
the 54-atom supercell. Again, this is in agreement with 
the "negligible" relaxation for the T site in silicon found 
by Luchsinger et al. 



B. Relaxations including zero— point motion 

The influence of the zero-point motion of the 
muon/proton on the relaxation of the silicon/germanium 
host lattice was studied by calculating the total energy, 
E"^, of Eq. 2.4, for different r elaxa tions of the NN host 



atoms, as described in Section [I C 



For the BC site four different relaxations of the NNs 
were considered, and for each of these the six NNNs were 
also relaxed. The NN relaxations are in addition to the 
static relaxations given in Section III A. The inclusion 



of the zero-point energy of the muon was found to give 
only a small correction to the static relaxations; the NN 
silicon atoms relaxed outwards by just an additional 0.01 
A in the [1 1 1] direction, so that the final separation 
of the muon from a NN silicon atom is 1.58 A in the 
16-atom supercell. The much smaller zero-point energy 
of the heavier hydrogen impurity is swamped by the in- 
crease in energy of the crystal as the separation of the NN 
atoms is increased and thus there is no additional relax- 
ation. As a check on the finite-size errors, the energies of 
five geometries were recalculated using the 54-atom su- 
percell. These energies and the corresponding potential 
energy curves calculated within the 16-atom supercell 
are shown in Fig. ^ The very small differences between 
the 16-atom and 54-atom results justify the use of the 
16-atom supercell in calculations of the shape of the po- 
tential well at the BC site. 

The story is very similar for the BC site in germa- 
nium. When the zero-point energy of the muon is in- 
cluded, the relaxation of the NN atoms again increases 
by just 0.01 A, so that the final separation of the muon 
from a NN germanium atom is 1.69 A in the 16-atom 
supercell. Once more, the smaller zero-point energy of 
the proton means there is no additional relaxation due to 
quantum effects. Fig. || shows the potential energy well 
and calculated wave functions for the muon and proton 
at the BC site in germanium. 



C. Zero— point energies 

The zero-point energy of the muon at the BC site was 
calculated to be 0.63 eV in silicon and 0.56 eV in ger- 
manium. It is perhaps surprising that such large zero- 
point energies have so little effect on the relaxations. As 
shown in Fig. ^, the potential well is narrow in the direc- 
tion along the bond and wider perpendicular to the bond. 
Within the harmonic approximation one can decompose 
the zero-point energy into contributions from the well 
along and perpendicular to the bond. For a muon at the 
BC site in silicon this gives 0.47 eV in the direction along 
the bond and 0.22 eV perpendicular to the bond. (The 
sum of these differs from the full zero-point energy of 
0.63 eV because the latter does not assume the harmonic 
approximation.) The corresponding energies for germa- 
nium are 0.37 eV and 0.22 eV in the directions along 
and perpendicular to the bond, respectively. If we were 
to consider only the zero-point energy in the direction 
along the bond then the outwards relaxation of the sil- 
icon/germanium atoms would be larger; approximately 
0.03 A in silicon and more than 0.025 A m germanium. 
Although the component of the zero-point energy along 
the bond is significantly reduced by further outward re- 
laxation of the silicon/germanium atoms, the potential 
well also gets narrower in the plane perpendicular to the 
bond, which tends to increase the zero-point energy. The 
narrowing of the potential well in the plane perpendicular 
to the bond correlates with the narrowing of the bonding 
charge cloud as the bond lengthens. 

Our result of 0.63 eV for the zero-point energy of the 
muon at the BC site in silicon is close to the value of 
0.54 eV obtained by Claxton et al. ||2^ from Hartree- 
Fock calculations on Si26ll3o clusters. In that calcula- 
tion the potential well at the BC site was assumed to be 
cylindrically symmetric about the bond (as it is in this 
work) and the resulting Schrodinger equation was solved 
within the harmonic approximation. 

The larger mass of the proton significantly reduces the 
quantum effects. We calculated the zero-point energy 
of a proton at the BC site to be 0.20 eV in silicon and 
0.18 eV in germanium. Our value for silicon is close to 
that of 0.18 eV obtained by Luchsinger et al. [Q, which 
suggests that the harmonic approximation to the poten- 
tial well used in that work is quite good for the proton 
ground state. 

In contrast to the results of plane-wave pseudopoten- 
tial calculations [0j28| , we find the T site corresponds to 
a local minimum in the potential energy surface. This 
is, however, in agreement with a more recent plane-wave 
pseudopotential study. The calculated energy surface 
along the [1 1 1] direction is shown in Fig. It turns 
out that the muon/proton is not strongly bound in our 
potential well, which turns over at the hexagonal site sit- 
uated at a distance of 1.18 A from the T site along the 
[1 1 1] direction. To confine the muon/proton in the well 
we therefore constrained the fit to prevent the potential 
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turning over, as shown in Fig. ^. The zero-point energy 
of the muon/proton calculated in such a well is then an 
upper bound on the true value, but as the wave function 
decays quite rapidly away from the T site this bound is 
accurate. 

In contrast to the BC site, investigation of the finite 
size effects present in the calculation of this energy sur- 
face showed that while the results in the 16- and 54-atom 
supercells were qualitatively similar, they differed signif- 
icantly in the openness and depth of the potential well. 
As a result, the calculation of the zero-point energy etc. 
of the muon/proton at this site was carried out using 
the potential well obtained from the 54-atom supercell. 
The expense of the LSDA calculations with this super- 
cell made the generation of data points off the [1 1 1] axis 
too costly. Therefore, the three parameters left undeter- 
mined after the one-dimensional fit to the data on the [1 
1 1] axis were assigned the values obtained in the fit to 
the 16-atom supercell data. This is not a critical choice 
since the one-dimensional fit has already constrained the 
shape of the energy surface in eight directions (due to the 
symmetry of the T site). This procedure was also used 
to generate a three-dimensional potential-energy func- 
tion from the fit to data along the [1 1 1] axis (calculated 
using the 16-atom supercell) at the T site in germanium. 
The parameters obtained from these fits are given in Ta- 
ble ^ 

For a muon at the T site the zero-point energy was 
calculated to be 0.28 eV in silicon and 0.22 eV in germa- 
nium. For a proton the corresponding values are 0.09 eV 
and 0.06 eV. The ground state wave functions of the 
muon/proton along the [1 1 1] axis in silicon are shown 
in Fig. 1^ and those in germanium in Fig. ^. The results 
for germanium must be considered approximate since we 
have not calculated any data points off the [1 1 1] axis 
in this case. In addition, the 16-atom supercell was used 
for all of the germanium calculations and therefore it fol- 
lows from the behaviour found in silicon that the true 
potential energy surface will be more open than the one 
we have obtained. 



bound within the well. In germanium, the potential well 
at the BC site is 1.51 eV deep. The first excited state is 
two-fold degenerate with an energy of 0.78 eV and is of 
the same character as in silicon. 

For the more massive proton the excited states are of 
lower energy. At the BC site in silicon, the two-fold 
degenerate first excited state of the proton has an energy 
of 0.27 eV, which is 0.07 eV higher than the ground state, 
while in germanium the excited state has an energy of 
0.25 eV, which is also 0.07 eV above the ground state. 

Each excited state of the muon/proton defines a dif- 
ferent adiabatic potential for the nuclei [i.e., a different 
£'^(r„) in Eq. 2.7). It is therefore possible for the lat- 
tice relaxations that occur when the muon/proton is in 
its first excited state (say) to be different from those for 
the ground state. For instance, the fact that the wave 
fimction of the first excited state of the muon/proton is 
essentially an excitation in the plane perpendicular to 
the bond, combined with the fact that the potential well 
in this plane becomes narrower as the separation of the 
NN atoms increases, results in the NN atoms actually re- 
laxing towards the impurity. This relaxation is small for 
the muon, but effectively zero for the proton due to the 
smaller zero-point energy. The effect on the energies of 
the excited states is negligible. 

At the T site in silicon, the potential well is consider- 
ably shallower with a depth of only 0.20 eV. For the muon 
this means that even the ground state energy (0.28 eV) of 
our constrained potential well (which is an upper bound 
on the true ground state energy) is greater than the well 
depth. The proton, however, has a (triply degenerate) 
first excited state with an energy of 0.14 eV which may 
therefore be bound at the T site. In germanium the po- 
tential well at the T site has a depth of just 0.18 eV 
in our 16-atom supercell calculations. It follows from 
the behaviour found in going from the 16-atom to the 
54-atom supercell in silicon that the true well depth in 
germanium will probably be less than this. It is there- 
fore unlikely that excited states of either the muon or the 
proton will be bound at the T site in germanium. 



D. Excited states of the muon and proton 

For a muon at the BC site in silicon our zero-point 
energy of 0.63 eV is considerably smaller than the well 
depth of 1.37 eV, indicating the possibility that excited 
states of the muon may be bound within the well. Numer- 
ical calculations show a two-fold degenerate first excited 
state at an energy of 0.84 eV. The wave functions of these 
states are similar to those obtained from a harmonic ap- 
proximation, i.e., they consist essentially of an excitation 
within the plane perpendicular to the bond. The energy 
of the first excited state is also reasonably well described 
within the harmonic approximation which predicts the 
excited state to be 0.22 eV above the ground state. It 
is also possible that some of the higher energy states are 



E. Energy barriers at the T and BC sites 

The heights of the energy barriers confining the muon 
and proton at the BC and T sites are clearly of great im- 
portance in determining the dynamics of the impurities 
within the lattice and hence are a significant part of the 
configuration-coordinate diagram. 

The static barriers (i.e., excluding zero-point effects) 
experienced by the muon and proton are identical. For 
the BC site in silicon we calculate the static barrier to 
motion towards the hexagonal site (in the [—1 1 0] di- 
rection) to be 1.37 eV while in germanium it is 1.51 eV. 
The effective barrier height is reduced by the zero-point 
energy and therefore depends on the nature of the im- 
purity. Including this effect, the effective barrier experi- 
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enced by a muon at the BC site in silicon is 0.74 eV while 
in germanium it is 0.95 eV. For the proton, the effective 
barriers (1.17 eV in silicon and 1.33 eV in germanium) 
are higher. The effective barrier height for the muon at 
the BC site may be considered a measure of the barrier 
to the BC^T site transition. In reality this transition is 
believed to involve charged states: muonium at the BC 
site is first ionized (with activation energy 0.22 eV |30|) 
and then moves to the T site while simultaneously recap- 
turing an electron to regain its neutral charge state. The 
sum of the activation and barrier energies for these two 
processes as measured experimentally is 0.60 eV. 

At the T site the energy barriers are very much lower. 
In silicon the static barrier to motion of the muon to- 
wards the hexagonal site (in the [1 1 1] direction) is cal- 
culated to be 0.20 eV while in germanium it is 0.18 eV. 
When zero-point effects are taken into account, the ef- 
fective barriers for the muon at the T site in silicon and 
germanium are zero indicating that, even at T = OK, 
the muon is free to diffuse through the interstitial region. 
However, this barrier is not appropriate for the T^BC 
site transition because in our calculations for the muon in 
the interstitial region, the host atoms around the BC site 
are unrelaxed. The process by which these atoms relax 
the large distances required to allow the muon to move 
to the BC site is unclear. Experimentally, the barrier for 
the T— >BC site transition in silicon is 0.39 eV. 

Since the zero-point energy of the proton is around 
a third of that of the muon, these calculations indicate 
that it will be bound at the T site in both silicon and 
germanium with an effective barrier of around 0.12 eV 
in each case. As previously discussed the true effective 
barrier in germanium will probably be lower than this. 



F. Hyperfine and superhyperfine parameters 

The hyperfine parameters depend on the spin density 
in the region at and around the atomic nuclei. More 
specific insight into the origin of the large measured dif- 
ferences between hyperfine parameters for niuons located 
at the two impurity sites can be gained from a consider- 
ation of the spin density isosurfaccs. Fig. ^ shows spin 
density contour plots for silicon in appropriate planes en- 
compassing the BC and T sites. Evidently the majority 
spin density around an impurity placed at the BC site 
is largely dispersed onto the two nearest-neighbour sili- 
con atoms; the spin density in a small region around the 
hydrogen nucleus is comparatively small and of opposite 
sign. At the T site, by contrast, almost all of the ma- 
jority spin density is localized on the defect. From these 
calculations one therefore expects the isotropic hyperfine 
parameter at the two sites to be of opposite sign, with the 
magnitude of the parameter at the BC site much smaller 
than at the T site. 

It is of course necessary to check the dependence of cal- 
culated hyperfine and superhyperfine parameters on the 



supercell size and basis set quality. A set of computed 
numbers are shown in Table III. The parameters ap- 



pear to be reasonably well converged with respect to the 
basis set, but the convergence with increasing supercell 
size is less good, particularly for the isotropic hyperfine 
and superhyperfine parameters at the BC site. Table IV 
gives the hyperfine and superhyperfine parameters cal- 
culated at the BC site in both silicon and germanium 
together with the results of other calculations for com- 
parison. Without motion averaging, the values obtained 
for silicon using the LSDA approximation are in reason- 
able agreement with both experiment and other DFT 
calculations. 

Both the hyperfine and superhyperfine motion- 
averaged tensors ((A^^) and (A^^^)) were found to be 
axially symmetric about the Si-Si bond ([1 1 1] direc- 
tion) in agreement with the experimental results. As 
Luchsinger et al. found, motion averaging increases 
the values of all but one (the anisotropic hyperfine pa- 
rameter) of the hyperfine and superhyperfine parameters, 
with the isotropic (contact) term on the muon being the 
most sensitive. This is because of the very small contact 
charge density which varies quite significantly with the 
muon position (Fig. H). In agreement with Luchsinger 
et al. 0, use of the Perdew-Wang |16| GGA functional 
did not consistently improve the values of the parame- 
ters. The results obtained for the muon at the BC site 
in germanium follow a similar pattern. 

The calculated hyperfine parameters for the T site are 
given in Table |v|. For silicon our values are in good agree- 
ment with both experiment and previous calculations. 
Again use of the Perdew-Wang ||l^ GGA functional fails 
to improve this agreement. For germanium, our value 
of the isotropic hyperfine parameter at the T site also 
agrees quite well with the measured value. 

The behaviour of the isotropic hyperfine parameter 
along the [1 1 1] axis in the vicinity of the T site in 
silicon and germanium is shown in Fig. Motion av- 
eraging for the muon/proton at the T site reduces the 
isotropic hyperfine parameter in both silicon and germa- 
nium. The final motion-averaged results are in reason- 
able agreement with experiment. Motion averaging of 
AJ^ resulted in an isotropic tensor, in agree ment with 



the symmetry arguments presented in Section II D and 
experimental observations. 

In a recent application of the path-integral Monte 
Carlo approach, Miyake et al. ||ll[] studied hydrogen and 
muonium at the T site in silicon, with the electron- 
electron interactions calculated within the LDA. They 
found the T site to be a local maximum in the po- 
tential energy surface, in agreement with Luchsinger et 
al. but in disagreement with our results and a recent 
plane-wave pseudopotential calculation. ]29| ] Their path- 
integral Monte Carlo study showed that quantum effects 
led to the muonium distribution being centred on the T 
site while hydrogen behaved as a largely classical particle 
and was thus distributed away from the local maximum 
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on that site. Evaluating the motion-averaged isotropic 
hyperfine parameter with our hyperfine data gives a value 
of 492 MHz for the hydrogen distribution of Miyake et 
at, but 685 MHz with our hydrogen distribution. There- 
fore if one could measure the isotropic hyperfine signal of 
hydrogen at the T site, one could deduce whether the T 
site is a maximum or minimum in the potential energy 
surface. 



G. Energies of a muon/proton at the T and BC sites 

The question of the relative stabilities of the muon and 
proton at the BC and T sites is of considerable interest. 
For a particular impurity this energy difference is the 
sum of contributions from the static-lattice energy and 
the zero-point energy. The contribution from the static 
lattice is sensitive to the size of the supercell and to the 
quality of the basis set. We investigated this point using 
the 16- and 54-atom supercells. We have added a BSSE 
correction to each of the static-lattice energy differences 
quoted here. With the 16-atom silicon cell and the stan- 
dard basis set, the T site was found to be 0.63 eV lower 
in energy than the BC site. Using the large basis set re- 
duced this to 0.32 eV. In the 54-atom supercell and using 
the standard basis set the T site was 0.41 eV lower in en- 
ergy than the BC site. With the large basis set, this was 
reduced to just 0.07 eV. These results indicate that a 16- 
atom supercell is too small to give reliable estimates of 
the static-lattice energy difference between the two sites. 
A summary of the computed energies that influence the 



relative stabilities is given in Table Vl 



There have been several previous calculations of the 
static-lattice energy difference between the T and BC 
sites in silicon. Using a plane-wave pseudopotential 
method and the LSDA, Chang and Chadi ||] found the 
T site to be lower in energy, but only by an amount 
<0.25 eV. Luchsinger et al. [Q, also using a plane-wave 
pseudopotential method, found the T site to be 0.15 eV 
higher in energy than the BC site within the LSDA and 
0.19 eV higher within the GGA. Note, however, that 
Luchsinger et al. 0| found the T site to be a local maxi- 
mum in the energy and that a nearby site has an energy 
about 0.05 eV lower. It is clear from the various results 
that the static-lattice energy difference between the T 
and BC sites in silicon is small within the LSDA/GGA, 
but its precise value has yet to be settled. 

The fact that the static-lattice energy difference is 
small means that the zero-point energy of the impurity 
is crucial in determining the relative stability of the T 
and BC sites. For a muon in silicon we have found the 
zero-point energy at the BC site to be 0.35 cV larger 
than at the T site. This difference is large enough to to 
make the BC site unfavourable for the muon, irrespective 
of which of the above values for the static-lattice energy 
difference is used. However, the zero-point energy of the 
proton at the BC site in silicon is only 0.12 eV higher 



than at the T site. Therefore, for this impurity the rela- 
tive stability of the two sites depends on the precise value 
of the static-lattice energy difference. 

In germanium with a 16-atom supercell, the difference 
in static lattice energies favours the T site by an energy of 
0.57 eV. The convergence with respect to supercell size 
found in silicon suggests that this difference in a fully 
converged LSDA calculation would be smaller. We es- 
timate that the zero-point energy of a muon at the BC 
site is 0.34 eV larger than at the T site (where the form 
of the potential energy surface was obtained from data 
calculated along the [1 1 1] axis only). For a proton the 
corresponding value is 0.12 eV. These results are similar 
to those obtained in silicon and thus it is likely that for 
the muon the T site is lower in energy. Without a fully 
converged value for the static lattice energy difference we 
are unable to draw any conclusions on the lowest energy 
site of the proton. 



IV. CONCLUSIONS 

We have calculated the zero-point motions and ener- 
gies as well as the hyperfine parameters of both muonium 
and hydrogen when present as impurities in silicon and 
germanium crystals at the BC and T sites. The electron, 
muon/proton and ion motions were decoupled using a 
double adiabatic approximation, and for the BC site we 
have included the effect of the zero-point motion on the 
relaxation of the host lattice. The ground states of both 
the muon and proton at the BC sites of silicon and ger- 
manium are strongly confined within a potential well of 
depth 1.37 eV (silicon) and 1.51 eV (germanium). The 
zero-point energy of a muon at the BC site is calculated 
to be 0.63 eV for silicon and 0.56 eV for germanium. De- 
spite the relatively large zero-point energy of the muon at 
the BC site, it causes only a small additional outwards 
relaxation of the nearest-neighbour silicon/germanium 
atoms of about 0.01 A. For the proton the additional 
relaxations of the nearest neighbours due to zero-point 
motion are negligible. At the T site the static relaxations 
of the host atoms are very small and the zero-point en- 
ergy is considerably smaller than at the BC site, being 
0.28 eV (0.22 eV) for a muon in silicon (germanium). 
It is therefore reasonable to assume that the additional 
relaxation due to the zero-point motion is negligible for 
either a muon or proton at the T site. 

The relaxation of the crystal around either the BC or 
T sites is practically independent of whether the impu- 
rity is a muon or a proton. This result confirms one of 
the underlying assumptions of the widely accepted con- 
figuration coordinate model. |^ The potential well at the 
BC sites of both silicon and germanium is reasonably well 
described by a harmonic approximation, at least for the 
ground states of the muon and proton. The potential well 
at the BC site in both materials is deep enough to bind 
several excited states of the muon and proton, although 
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we are not aware of any experimental evidence for such 
states. The potential well at the T site in either silicon or 
germanium is not deep enough to bind the muon which 
is free to diffuse through the interstitial region, although 
our calculations suggest that the proton is bound at this 
site at T = K. 

Various LSDA and GGA calculations have indicated 
that the energies for a static muon or proton at the BC 
and T sites in silicon are very similar. However, we 
have calculated the difference in zero-point energies of 
a muon at the T and BC sites in silicon (germanium) to 
be 0.35 eV (0.34 eV) which is sufficient to make the T site 
more stable, whether we assume our value for the static- 
lattice energy difference between the BC and T sites or 
those of others. [^,0 This result is in conflict with the 
interpretation of experimental data. 

The hyperfine parameters calculated for silicon in our 
all-electron calculations are close to those obtained in 
plane- wave pseudopotential calculations. This agree- 
ment confirms that the procedure used to correct for the 
pseudopotential and for the incomplete plane-wave basis 
sets are accurate. For silicon our static LSDA results are 
in reasonable agreement with other LSDA results and ex- 
periment. Our hyperfine parameter for the muon at the 
T site in germanium is in much better agreement with 
experiment than the only previous calculation of which 
we are aware. |^ 

In our work the motion averages of the hyperfine and 
superhyperfine parameters are evaluated by averaging 
over the squared modulus of the wave function obtained 
from the full solution of the muon/proton Schrodinger 
equation in the potential well. We note that the symme- 
try of the potential wells requires that the symmetry of 
the motion-averaged hyperfine tensors at the T and BC 
sites are the same as if the muon/proton was situated ex- 
actly at the sites. We have obtained detailed information 
about the variation of the hyperfine and superhyperfine 
parameters with the position of the muon/proton. Our 
results show that motion averaging for the muon/proton 
at the BC site in silicon and germanium increases the 
values of all of the hyperfine and superhyperfine param- 
eters apart from the anisotropic hyperfine term which 
decreases slightly, in agreement with the conclusions of 
Luchsinger et al. Q With the exception of the isotropic 
hyperfine term however, all of the changes are small. At 
the T sites in silicon and germanium, motion averaging 
reduces the isotropic hyperfine parameter. 
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TABLE I. The values of the parameters in Eq. 2.8 defining 
the potential well at the BC site of silicon and germanium. 
The fit is applicable within a cylinder centred on the BC site 
of radius 1.0 A (1.1 A) in the plane perpendicular to the bond 
and up to ±0.5 A (0.37 A) along the direction of the bond in 
silicon (germanium). The units are such that if the lengths in 



Eq. I.i are expressed in Bohr radii then the potential energy 
is in Ha. N.B. The values quoted are those carried into the 
zero-point calculation and therefore the number of significant 
figures should not be taken as an indication of the accuracy 
of the fit. 

Silicon Germanium 
0.00588753 0.00610711 
7 0.0807689 0.0548517 
<5 0.0 0.0 
C 0.000337036 0.000224779 
r) 0.0654643 0.0506238 



T ABL E IL The parameters defining the expansions 
(Eq. |2.9| ) of the potential energy surface around the T site 
in silicon and germanium. The fit is applicable over the re- 
gion bounded by a sphere of radius 1.0 A centred on the T 
site. The units are such that if the lengths in Eq. 2.£ are ex- 
pressed in Bohr radii then the potential energy is in Ha. N.B. 
The values quoted are those carried into the zero-point calcu- 
lation and therefore the number of significant figures should 
not be taken as an indication of the accuracy of the fit. 



Silicon 



Germanium 



a2 
as, 

a4, 

&4 

as 
ae 
h 

C6 



0.00440492 

0.00589898 
-0.00197527 

0.000930883 
-0.00456177 

0.00155753 
-0.000556122 

0.00696501 



0.00151269 
0.00542489 
0.000339949 
-0.000160207 
-0.00311144 
0.000775975 
-0.000277064 
0.00347002 
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TABLE III. The dependence of hyperfine and superhyperfine parameters for muonium in silicon 
on supercell size and basis set. 



BC site T site 



Supercell 


Basis set 


As, 








As, 


16 atom 


Standard 


-27.1 


17.7 


-147 


-13.6 


2302 


16 atom 


Large 


-21.4 


13.0 


-114 


-10.7 


2366 


54 atom 


Standard 


-1.6 


16.4 


-91.0 


-12.4 


2362 


54 atom 


Large 


4.5 


9.8 


-57.1 


-8.1 


2389 



TABLE IV. Static and motion-averaged (indicated by {)) hyperfine parameters for the muon at 
the BC site and the nearest-neighbour atoms. PS denotes a pseudopotential calculation. 



Hyperfine parameters (MHz) 
Silicon Germanium 





As, 


Ap, 


As^i 




As, 


Ap, 


Asc„ 


A 

^PGe 


LSDA'' 


-1.6 


16.4 


-91.0 


-12.4 










LSDA'' 


-27.1 


17.7 


-147 


-13.6 


-24.6 


16.4 


-80.7 


-5.8 


(LSDA)'' 
GGA*' 


2.5 


14.5 


-141 


-13.4 


3.6 


12.5 


-75.5 


-5.6 


-89.3 


18.9 


-155 


-14.0 


-64.6 


17.1 


-81.0 


-5.9 


LSDA'^ 


-104 


58.5 


-127 


-53.5 


-87 


64 


-85 


-24 


PS-LSDA'^ 


-26 


22.8 


-90 


-20.2 










PS-GGA'^ 


-81 


27.5 


-192 


-28 










(PS-GGA)'^ 


-65 


21.7 


-191 


-26.2 










PS-LSDA^ 


-26.8 


18.1 


-83.8 


-22.7 










PS-LSDA' 


-35 


22.3 


-85 


-21.5 










Experiment^ 


-67.3 


25.3 


-95.1 


-21.2 











Experiment'' -96.5 34.6 



^This work, 54-atom supercell and "standard" basis set. 

^This work, 16-atom supercell and "standard" basis set. 

'^Casarin et al. p7|] 

'^Luchsinger et al. H] 

°Van de Walle. g 

fVan de Walle and Blochl. @ 

SKiefl and Estle. |^ 

'^Patterson. § 
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TABLE V. Isotropic hyperfine parameters for the muon at the T site of silicon and germanium. 
The quoted resuhs are for the "standard" basis set. 





bihcon 


As, (MHz) 


Germanium 


LoUA 


2302 




2236 




2651 




2548 


(LSDA)" 
LSDA'^ 


2096 
2362 




2032 


(LSDA)'' 
PS-LSDA'^ 


2152 
1939 






PS-GGA" 


2098 






LSD-VBH'^ 


3043 




3977 


Experiment"^ 


2006 




2360 



^This work, 16-atom supercell. 
^This work, 54-atom supercell. 
"^Luchsinger et al. 
'^Casarin et al. [|^, extended basis set. 
''Patterson ^ and references therein. 



TABLE VL A summary of the energies influencing the relative stability of the BC and T sites 
in silicon and germanium. The importance of the zero-point energy of the muon in determining 
the favoured site is clear. 

Silicon Germanium 





BC site 


T site 


BC site 


T site 


Static lattice energy w.r.t. BC site (eV) 





-0.07=^ 





-0.57" 


Muon zero-point energy (eV) 


o.es"^ 


0.28'' 


0.56*^ 


0.22'' 


Total energy w.r.t. BC site for muon (eV) 





-0.42 





-0.91 


Proton zero-point energy (eV) 


0.20"^ 


cog"* 


0.18*^ 


0.06'' 


Total energy w.r.t. BC site for proton (eV) 





-0.18 





-0.69 



^54-atom supercell and "large" basis set. 
''16-atom supercell. 

''16-atom supercell and "standard" basis set. 
''54-atom supercell and "standard" basis set. 
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FIG. 1. The muon/proton at the bond-centred (a) and tetrahedral (b) sites in silicon. The 
dashed circles in (a) show the unrelaxed positions of the silicon atoms. These are not shown in (b) 
because the relaxation around the T site is negligible. 




FIG. 2. The square of the proton (thin solid line) and muon (dashed line) wave functions at the 
BC site in silicon. The symbols denote the potential well for the 16- and 54-atom supercells, and 
the thick solid line is the fit of the 16-atom data to Eq. 2.8 with the parameters of Table |. 
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FIG. 3. The square of the proton (thin sohd hne) and muon (dashed line) wave functions at the 
BC site in germanium. The sy mbo ls denote the potential well for the 16-atom supercell, and the 
thick solid line is the fit to Eq. 2.8 with the parameters of Table \A 




FIG. 4. The calculated potential experienced by the muon/proton in the {110} planes for the 
BC site in silicon. The figure shows both of the NN and two of the NNN silicon atoms of the 
muon/proton with the bond lengths drawn to scale. The contours range from 0.1 to 0.7 eV in 
increments of 0.1 eV. 
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FIG. 5. The square of the proton (thin sohd hne) and muon (dashed line) wave functions at the 
T site in silicon. The symbols denote the potential well for the 54-atom supercell, and the thick 
solid line is the fit to Eq. 2.£ with the parameters of Table The potential has a maximum at the 
hexagonal site situated at 1.18 A from the T site, but the fit has been constrained so that it forms 
a simple potential well. 
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FIG. 6. The square of the proton (thin sohd hne) and muon (dashed line) wave functions at the 
T site in germanium. The sy mbo ls denote the potential well for the 16-atom supercell, and the 
thick solid line is the fit to Eq. 2.9 with the parameters of Table ^ The potential has a maximum 
at the hexagonal site situated at 1.23 A from the T site, but the fit has been constrained so that it 
forms a simple potential well. 
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FIG. 7. Spin density contour map in the neighbourhood of (a) the BC site and (b) the T site in 
sihcon. In (a) the muon position is located dead centre, with the two nearest-neighbour sihcon atoms 
above and below. In (b) the large concentration of positive spin density is located on the muon 
position, the positions of nearby silicon atoms in this plane are indicated with crosses. Continuous, 
dashed and dot-dashed lines correspond to positive, negative, and zero values respectively. The 
separation between adjacent isodensity contours is 0.001 e/bohr^ 



(8) 








X 











(a) (b) 



FIG. 8. Variation of the isotropic hyperfine parameter and the xy component of Ap with dis- 
placement of the muon/proton from the BC site. The calculations were performed in silicon with 
the 16-atom supercell and the standard basis set. The lines are guides to the eye. 
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FIG. 9. The variation in the isotropic hyperfine parameter along the [1 1 1] direction at the T 
site in sihcon (54-atom supercell) and germanium (16-atom supercell). A positive displacement 
indicates movement towards the hexagonal site. The lines are guides to the eye. 
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